### This script will intersect BRT predictions with the dataframe used to build the model
### The resulting data frame will contain for all site/days Obs, Preds, Resids, and Topographic Data used to construct the model

### Establish directories

minpreds.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/brtpreds/tmingzip/'
maxpreds.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/brtpreds/tmaxgzip/'
out.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/'

### Load library

library('SDMTools')

### Read in .csv file with data used to build the model
### This will be the framework for the resultant dataframe

model.data = read.csv('/home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/final_new.csv',header=T)

### Remove rows from model.data where micro_min or micro_max is 'NA'

model.data = model.data[which(is.na(model.data$micro_min)==F),]

### Get a list of the ASCII files that are in the preds directory

minpreds.list = list.files(minpreds.dir,pattern='.asc.gz',full.names=T,recursive=T)
maxpreds.list = list.files(maxpreds.dir,pattern='.asc.gz',full.names=T,recursive=T)

### Use a for loop

### Blank object to bind data onto

out.data = NULL 

### Unique dates to cycle through

date.list = c(unique(model.data$date))

### Position tracker

i=1

for (datex in date.list)

	{
	
	# Read in the microclimate predictions for the day corresponding to datex
	
	max.asc = read.asc.gz(maxpreds.list[which(substr(maxpreds.list,92,99)==datex)])
	min.asc = read.asc.gz(minpreds.list[which(substr(minpreds.list,92,99)==datex)])
	
	# Subset model.data by datex
	
	t.data = model.data[which(model.data$date==datex),]
	
	# Extract preds for all sites in t.data
	
	t.preds = data.frame(t.data,maxpreds=extract.data(cbind(t.data$lat,t.data$long),max.asc),minpreds=extract.data(cbind(t.data$lat,t.data$long),min.asc))
	
	# Re-assemble into one dataframe for writing out
	
	out.data = rbind(out.data,t.preds)
	
	cat(i/length(date.list),'\n')
	
	i=i+1
	
	}
	
# Close loop

write.csv(out.data,file=paste(out.dir,'BRT_Model_Obs_&_Preds.csv',sep=''),row.names=F)

# Done
	
	